{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 从疝气病症预测病马的死亡率"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "导入模块"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import *\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "载入数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loadDatSet(filename):\n",
    "    dataMat = [] # 用于存放数据\n",
    "    labelMat = [] # 用于存放标签\n",
    "    fr = open(filename) # 打开文件\n",
    "    for line in fr.readlines(): # 读取行\n",
    "        lineArr = line.strip().split() # 拆分\n",
    "        dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])]) # 整理数据\n",
    "        labelMat.append(int(lineArr[2])) # 加入标签\n",
    "    return dataMat,labelMat # 返回数据和标签"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "sigmoid函数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(inX):\n",
    "    return 1.0/(1.0+exp(-inX))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "计算梯度"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calcGrad(weights,dataMatrix,labelMatrix):\n",
    "    h = sigmoid(dataMatrix*weights) # 计算sigmoid输出\n",
    "    grad = -dataMatrix.T*(labelMatrix - h) # 求梯度\n",
    "    return grad # 返回梯度"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "梯度下降更新算法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradientDescent(dataMat,labelMat):\n",
    "    dataMatrix = mat(dataMat) # 将数据转换成矩阵形式\n",
    "    labelMatrix = mat(labelMat).T # 将标签转换成列矩阵\n",
    "    n = shape(dataMatrix) # 得到数据矩阵的尺寸\n",
    "    weights = ones((n[1],1)) # 初始化回归系数\n",
    "    alpha = 0.001 # 设置学习率\n",
    "    echo = 500 # 设置迭代次数\n",
    "    for i in range(echo):\n",
    "        weights = weights - alpha*calcGrad(weights,dataMatrix,labelMatrix) # 更新回归系数\n",
    "    return weights # 返回回归系数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "随机梯度下降"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stocGradDescent(dataMat,labelMat,numIter=150):\n",
    "    dataMatrix = mat(dataMat) # 将数据转换成矩阵形式\n",
    "    labelMatrix = mat(labelMat).T # 将标签转换成列矩阵\n",
    "    n = shape(dataMat) # 得到数据矩阵的尺寸\n",
    "    weights = ones((n[1],1)) # 初始化回归系数\n",
    "    for j in range(numIter):\n",
    "        dataIndex = list(range(n[0]))\n",
    "        for i in range(n[0]):\n",
    "            alpha = 4/(1.0 + j + i) + 0.01 # 设置学习率更新方式\n",
    "            randIndex = int(random.choice(dataIndex)) # 随机选取数据索引\n",
    "            weights = weights - alpha*calcGrad(weights,dataMatrix[randIndex],labelMatrix[randIndex]) # 更新回归系数\n",
    "            dataIndex.remove(randIndex) # 删除当前索引\n",
    "    return weights # 返回回归系数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "画拟合直线"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotBestFit(wei,dataMat,labelMat):\n",
    "    weights = array(wei)\n",
    "    dataArr = array(dataMat)\n",
    "    n = shape(dataArr)[0] # 得到数据量\n",
    "    x1cord1 = [];x2cord1 = [] # 存放1类横纵坐标\n",
    "    x1cord2 = [];x2cord2 = [] # 存放0类横纵坐标\n",
    "    # 存放坐标\n",
    "    for i in range(n):\n",
    "        if int(labelMat[i]) == 1:\n",
    "            x1cord1.append(dataArr[i,1])\n",
    "            x2cord1.append(dataArr[i,2])\n",
    "        else:\n",
    "            x1cord2.append(dataArr[i,1])\n",
    "            x2cord2.append(dataArr[i,2])\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111)\n",
    "    ax.scatter(x1cord1,x2cord1,s=30,c='green',marker='s')\n",
    "    ax.scatter(x1cord2,x2cord2,s=30,c='blue')\n",
    "    x = arange(-3.0,3.0,0.1) # 直线横坐标\n",
    "    y = (-weights[0]-weights[1]*x)/weights[2] # 直线纵坐标\n",
    "    ax.plot(x,y,c='red')\n",
    "    plt.xlabel('X1')\n",
    "    plt.ylabel('X2')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XuUXGWZ7/HvkwsgkGAgAQLpEJRIuOTC0Caog+CAmYAXGEcOoAyoaIAlHJlzPODlKIrOQlFHXTJcgjCCw+g5cwRlQVCizqzADAQCdieQcAkI6TYhNJdAQlDSyXP+2NVWdXVVde3u2nu/e9fvs1av7tq1q/JUV2c/9b7v876vuTsiIiLNGpN1ACIiki9KHCIiEosSh4iIxKLEISIisShxiIhILEocIiISixKHiIjEosQhIiKxKHGIiEgs47IOIAmTJ0/2GTNmZB2GiEhuPPTQQy+4+5Rmzi1k4pgxYwYrV67MOgwRkdwws2ebPVddVSIiEosSh4iIxKLEISIisSSeOMzsRjN73sweqTj2FTP7g5l1lb5OrvPYRWb2uJmtM7PPJR2riIgML40Wx4+ARTWOf9fd55W+llbfaWZjgX8CTgIOB840s8MTjVRERIaVeOJw9+XASyN46Hxgnbs/7e5vAD8FTmlpcCIiEluWYxwXmtmqUlfWpBr3Hwj0VNzuLR0TEZEMZZU4rgHeCswDNgLfqXGO1ThWd59bM1tsZivNbGVfX19rohTJsZ4euOgimD8/+t7TM/xjRJqRyQRAd9808LOZXQ/cUeO0XqCj4vY0YEOD51wCLAHo7OzURurS1np6YO5c2LoVtm+Hri645Rbo7oaOjuEfL9JIJi0OM5tacfNvgEdqnPYgMNPMDjazXYAzgNvTiE8k7668spw0IPq+dWt0XGS0Em9xmNlPgOOByWbWC1wGHG9m84i6np4BziudewDwQ3c/2d37zexC4FfAWOBGd3806XhFimDFinLSGLB9OzzwQDbxSLEknjjc/cwah2+oc+4G4OSK20uBIaW6ItLYggVR91Rl8hg/PhrvEBktzRwXKaBLLoE994ySBUTf99wzOi4yWkocIgXU0RENhJ93XtTKOO88DYxL6xRyWXWRoujpiQa0V6yIup8uuaT5i39HB/zgB8nGJ+1JiUMkUCqplVCpq0okUCqplVApcYgESiW1EiolDpGAVC4T8sc/wriqzuRQSmq1nEl70xiHSCCqxzTGjYMdO6Lv/f3hlNRq7EXU4hAJRPWYRn9/lDQOOyysklqNvYhaHCKBqDem8aY3RfeFQmMvohaHSCAWLCjP9B4QyphGpbzEKclR4hAJRF6WCclLnJIcJQ6RQORlmZA04lTVVtjMvXh7HnV2dvrKlSuzDkNERqC6amugRRNiEi0SM3vI3TubOVctDhEJiqq2wqfEISJBUdVW+JQ4RCQoqtoKnxKHiNSU1QB1klVbGnRvDQ2Oi8gQWQ9QD+xDcs89sHMnjBkDxx4bbz+SWs+pQff6ghocN7Mbzex5M3uk4ti3zOwxM1tlZreZ2ZvrPPYZM1ttZl1mpkwgkpKsB6g7OqIksX49PPZYdHG/7rrowj/SVkLWr6lI0uiq+hGwqOrYMuBId58DPAF8vsHj3+Pu85rNhCIy2Ei6Z0IYoG71hT6E11QUiScOd18OvFR17G537y/dvB+YlnQcIu1ooHvmuuvgwQeb/9QewgB1qy/0IbymoghhcPwTwF117nPgbjN7yMwWN3oSM1tsZivNbGVfX1/LgxTJoy99CTZvjv+pPYRlRVp9oQ/hNRVFponDzL4I9AO31DnlXe7+F8BJwKfN7N31nsvdl7h7p7t3TpkyJYFoJWmqeGmtnh74l3+B6vqXZj61h7D8Sasv9CG8pqJIparKzGYAd7j7kRXHzgHOB05w921NPMdXgK3u/u3hzlVVVf6o4qX1LroIrr46qkqqZAZnnw0TJkTdQQsWjK5aKUkD1VUPPBBd7EONswjiVFVlsh+HmS0CLgWOq5c0zGwPYIy7byn9vBC4PMUwJUWNBkJ/8INsY8urFSuGJg2IEsfPfw7btoW/g19Hh97/EKVRjvsT4D7gUDPrNbNzgauACcCyUqnttaVzDzCzpaWH7gfca2bdwAPAne7+y6TjlWyo4qX1ao0RjBkDM2eWkwaoLFXiS7zF4e5n1jh8Q51zNwAnl35+GpibYGgSkAULok++lckj5IqXgS6UkLt6LrkkaklUd//ttpuStIxOCFVVIrmqeBlpiWva6g0GH3usylJldLTkiAQjLwOhF10UJYvq1tF55+WjPz6kQoTqlttZZ0WVYCG35IoqzuC4Eoe0nZ6eaH7DXaXZQyedBF/7WvMXqPnzo5ZGreMrVrQuziSFkKSrE9i4cbBjB4wdC/39qqxLW/BVVSJZ6emB2bPhlVfKx266KaoyWr26uQtU3sZjagmhWqm6kq6/f/B3VdaFS2McUgjNTh688kp49dWhx7dsab6qKE/jMSGrVUlXTYP2YVKLQ3Kvusuj0byEFSuGzqSGaL5DsxeogUHnrLt68q5Wy61a3lpy7UItDsm9OKuoLlgQTYCrNmZMvAvUQFfPihXRdyWN+KpbbuPGRe/NuNLHWbXkwqXEIbkXZ/LgJZfAxIlDj0+Y0B4XqJDWA6suFz7/fLjvvui71pIKm7qqJPfiDFZ3dESD4I2qqkKe3Dea2Kq79H73O7j+epg1a/S7641UrUH6BQvSjUHiUzmu5F4r5yWENMeh1bHVmn8yIKTXKdkIautYkaS1crnskLcXHW1sjaqYQnqdEj4lDimEVg1WJ7HYYqvGFUYbW61FD+M8V0jjI5ItjXGIVGj15L44pcJJx1a96GG1Rs/Vytch+acWh+Raqz8Ft3pyXyu7vkYbW2WX3ty5sOuuzZe+htyFJ+lT4pDcSmKV2lZvLzqa7qXqpAijj22gS6+rC558svnSV+2XIpXUVSW5ldSugY3WcYpbDnv44fDQQ4N34mume6lR11Cr1m2Ks15VEdbnktZRi0NyK+1PwXFbOD090eKJ1du37r778N1LoXUNFWV9Lg3wt4YSh+RWrSqhJD8Fx72YX3lltEVrJTM49dThu5dC6xpqdRdeFvKyAVcepJI4zOxGM3vezB6pOLa3mS0zsydL3yfVeew5pXOeNLNz0ohX8iHtT8FxL+a1zneHtWuH/7fSTorNyPv6XKG14vIsrRbHj4BFVcc+B/zG3WcCvyndHsTM9gYuAxYA84HL6iUYaT9pfwqOezEfzcW/KF1DIQmtFZdnqSQOd18OvFR1+BTgptLPNwGn1njoXwPL3P0ld38ZWMbQBCRtLM1PwXEv5qO5+KeVFNupzz/EVlxepbZWlZnNAO5w9yNLtze7+5sr7n/Z3SdVPeazwG7u/vXS7S8Br7v7txv9W1qrSpISd8vVELZorSfkdbmS0G6vN64ibR1bY+cEamY6M1sMLAaYPn16kjFJG4u75WoaW7SOdMXcpMqZQ6UNuFony8SxycymuvtGM5sKPF/jnF7g+Irb04D/qPVk7r4EWAJRi6O1oYqEaTRLgYTQ599M0mvlMvch7LVeBFmW494ODFRJnQP8osY5vwIWmtmk0qD4wtIxCUw79ZWHpF6r4X3vG/69qLfo4euvp/P+NVMeqxLaQLl74l/AT4CNwHaiVsS5wD5E1VRPlr7vXTq3E/hhxWM/AawrfX28mX/v6KOPdknP+vXukya5jx/vDtH3SZOi45Kst789+p3X+2r0Xgy8b+PGDX7MuHHpvH8XXlj+m6mM98IL450jrQGs9Cav6WlVVZ3p7lPdfby7T3P3G9z9RXc/wd1nlr6/VDp3pbt/suKxN7r7IaWvf04jXolH9fHZaWap9FrvxUD3T0fH0K10+/vTef+a6Spbvjz77jQZKvTBccmBEPrK29VwS6XD0PeielykmcckYbj1r3p64PHHhz5u3DiV0GZNS47IqKk+PjvV8z1mzy4vlT6g+r2obiHWksb7N9w8lyuvhB07hj5u7FhNhMyaEoeMWr0LwFln5XvAPIQB/2ZiqJwEeeedMGFC40mHjbaQrfeYJAw3yXHFiqjbrNqsWSqhzVyzgyF5+tLgePrWr48GLOfPj77ff3++B8xDGPAfaQzV70X1+fUGnGfPrv+YLGhgPF3EGBzP/CKfxJcSR/by/p8+hPiTiiGEpNiMvMRZFHESh7qqJBF5HzAPIf6kYsjLEul5iTNV27ZFfwBLlkR/IBlRVZUkIu87xoUQf7MxjGRmdV5mUOclzkQ891yUKbu6yl9PPFHeGeySS6I3PAOpLXKYJi1ymL28LygXQvzNxBBCnDJKO3ZECaGra3Ci2LSpfM706TBv3uCvGTOincFapEiLHEpO5X1BuRDibyaGdluoMPe2bIFVq8oJorsbVq+O1nmBKPMfcQScdFKUHObOjb4mhbUNkVocIjk2f360hlOt4xl2gYs79PYObkF0d8O6deVz9t67nByOOir6PmsW7LJLrH9q4hUT2fLGliHHJ+wygVc//2rTz6MWh0ibaMVYTKsuPG1r+/ZoP+Dq8YiXKvauO+SQKDGcfXY5SUyb1pKuplrvXaPjraDEIZJj1UuOjGTyXhYXntx6+eXB3UxdXbBmDbzxRnT/brvBnDnwt39bbk3MmRPNyiwQJQ6RHAthLKaQ3OGZZwZ3M3V1wbPPls/Zd9+o9bBwYTlJvO1tQ9d8KaDiv0KRDLRy86HhtHXJaiv88Y/w6KODE0R3N7xa6qYbMyZKCO98J1xwQTlJ7L9/tnFnSIlDpMVGsyufJKyvb+hYxGOPlVdT3HPPqGvprLPKZa9HHAG7755t3IFR4hBpMZXIBmDHDnjqqaFdTRs2lM+ZNi1KDKeeWk4Sb3lL1MLIkQm7TKhb3JAUJQ6RFgthuZI4srjwtNRrr0VzISqTxKpV0fIcEI05HH44nHhi1BQc6GraZ59s426RLCrflDhEWiyE5UriyE3JrTts3Dh0hvWTT0b3Aey1V5QYPvWpPyeIyXe8mxd3rgJWwRbgnuhL5cYjl1niMLNDgf9TcegtwJfd/XsV5xwP/AL4fenQre5+eWpBioxAK0pki2BUBQLbt0fb/1UniRdeKJ9z8MFRcvjIR8qtiIMOGjI34sXbt9b8J1RuPHKZJQ53fxyYB2BmY4E/ALfVOPUed39/mrGJjIZKZGMWCLzyStS1VNnV9Mgj8Kc/Rffvums0QP3BD5bHIubMiVoXkolQuqpOAJ5y92eHPVMkB9q9RLZmgcAW5/r/vZ7LP1Q1YP3735cfOHlylHEuuqg8HnHooUP3JpZMhZI4zgB+Uue+d5hZN7AB+Ky7P5peWCIyEg/f9yeO2L6GuXQzjy7m0cXc/m4m3bwZbibqTpo5E97+dvjkJ8stialTW7riqyQj88RhZrsAHwQ+X+Puh4GD3H2rmZ0M/ByYWed5FgOLAaZPn55QtCIyxIsvDlmGY/kjaxhLtGH4a+zOKubwb2NOZ69j53L6N46C2bNhjz0yDlxGKvPEAZwEPOzum6rvcPdXK35eamZXm9lkd3+hxrlLgCUQrY6bZMCSnTRnZEuVnTvh6aeHzo3o7S2fc8ABMHcurx37Pv7+R3N54E9zWdM/k7Hjx0b7hPwYSPn9yn25cYBCSBxnUqebysz2Bza5u5vZfGAM8GKawUk4NCN75GKvgPv669EAdeUM61Wrol8+wNix0djDcccN3jdi332jfw/4yqVRkt894wIBldy2XqaJw8x2B94LnFdx7HwAd78W+DBwgZn1A68DZ3gRNxCRpmhG9sg1XAG3covSge+PP17eonTChCgpfOxjfx6LmHLn8bzga4A18Not8F/Afw1ORO1eIFBkmSYOd98G7FN17NqKn68Crko7rqIo2j4LeZuRnYTRvKdjdsLbXoR5z8Hc56Lv854DvjK1fNJBB0VJ4rTTyi2JGTOGLMPxwh2aG9HOQuiqkoQUbZ+FvM3ITkLT7+mWLYOW4bh/KczeBLtH49W8MQYe3RfumgkfP+e75SQR2BalEiYlDskNzciuwWHaq1ELgq9/ve4WpVv3gms7oXt/6Nof1k6G7aX//R+/+OJMQpf8UuKQ3Gj7GdnbtzO71L00d1O5q2mf1wdO+FJ5i9Jzzim3IqZN48TL87Xiq4RNiUNyJa0B18zLfge2KK0ctH70UVaVdih9fRw8si/celjUgujaH/7zu6/AxIk1n04lqdkr0pijEodIlVTLfqu3KB1IEpVblO63XxTQxRdz5lNX0r0fPLEP7Bhb9Vx1kga0viRViSi+Io05KnEUmP5zj0xiZb8DW5RWrvZaa4vSY46pu0XpnVdcE8R7mrdPyNJaShwFpv/cI9OSst++vqFLglduUbrHHlFS+OhHo+9HHQVHHjnsFqV6TyUEShwiVWKV/cbdonQgSQyzRWkR+sOL8BqkNiUOCUbmA9Il9cp+L73wNbh/9eAEsXp1tHUplLcoPeGEwctwTJ4cO4Yi9IcX4TVIbQ0Th5lNBKa4+1NVx+e4+6pEI5O2EtI6VB3TnNW/2sjtl3fhXd0cs1sXs3d2Mf6wqi1K584dvCT4YYdFmw5lQJ/uw1ekMce6icPM/hvwPeB5MxsPfMzdHyzd/SPgL5IPT9pFZutQ9ffX3KL0wL4+Lhg45+CDoyRx9kfKSWL69KD2jdCn+/AVKYE3anF8ATja3TeWVqb9sZl9wd1vBcL5HyOFkMo6VJVblA4kieotSo88MtqidGD3uTpblOoTfmMTr6hfGiz51yhxjHP3jQDu/oCZvQe4w8ymAVqhVlqqpetQucP69YMrmmptUTpv3oi3KNUn/Mb0eyi2RonjVTN768D4RqnlcTzRLnxHpBGctI8Rr0P1xhuwZs3Q0tfNm6P7q7coHahqCnyL0iL1h1crwmtod40Sx6VUdUm5+xYzW0TtbV5FRqypdahqbFHKmjXROAVEcyDmzIHTTy+PReR0i9Iid3el9drUnZicRonjJuA6M/uOu/cDmNl+wHeAQ4HLU4hP2sif16HauTPqVnqgC5ZUzI/o6SmfPHVq1Hp43/vKXU2HHBLtTNeGitxCGSl1JyanUeI4GrgC+J2ZfQaYDfwP4Erg7BRik3awbVs0QF3ZkujuHrxF6axZcOyx5VZExRalEtEnaElT3cTh7i8D55eSxq+BDcAx7t5b7zEiDW3aNHSGdb0tSgfGIo44AnbbLdOwa9En/Mb0+ym2RvM43gx8E1gALAJOBu4ys8+4+29Tik/yaMcOeOKJoeMRzz1XPmf69Kj1cNpp5a6mGluUhkqf8BvT76fYGnVVPQxcDXy6NMZxt5nNA642s2fd/cxWBGBmzwBbgB1Av7t3Vt1vwPeJEtc2oomID7fi35YWqNqilO7u6Pbrpd2Fxo+PWg2LFg1ehkNblMoIpDHgrUH14TVKHO+u7pZy9y7gnWb2qRbH8R53f6HOfScBM0tfC4BrSt8lTe7Q2zu0FVG5RemkSVFyOP/88njErFmwyy7ZxS1BGunFOc6A90i7yzSoPrxGYxx1xzLc/fpkwqnpFOBmd3fgfjN7s5lNHZicKAnYvh3Wrh06N+Kll8rnVG5ROtDVNG1a0HMjJBxpXJzVOkhOCKvjOlE3mAPXufuSqvsPBCrqMOktHRuUOMxsMbAYYPr06clFWzSbNw+dYb1mTTSxDqKB6dmz4UMfKrci5syJBrKlIXV5SFGFkDje5e4bzGxfYJmZPebuyyvur/URdsiSJ6WEswSgs7NTS6JUc4/mRlR3NVVuUbrvvlFiWLiw3Ip429ui5cIzkueLr7o8pKgyTxzuvqH0/Xkzuw2YD1Qmjl6gcv7wNKLSYKlnYIvSygRRa4vSd7wj2qJ0IElUbFEaCl18RcKTaeIwsz2AMaWlTPYAFjJ0RvrtwIVm9lOiQfFXNL5Roa9vaCti7draW5QOdDU1sUWpSIjSmB+iOSjDy7rFsR9wW1RxyzjgX939l2Z2PoC7XwssJSrFXUdUjvvxjGLNVvUWpQOJonKL0gMPjCbNnXJKufT1rW/NzdyIrOS5OyyvRnpxTuP90Hs+vEwTh7s/Dcytcfzaip8d+HSacWXutdeiuRCVg9bVW5Qedhj81V9FiWIgSeyzT7Zx55S6w9Kni3O+Zd3ikM2b4b77Bnc1PfHE4C1K582Dc88tdzUdfnhmW5RK89TlIUWlxJG1e++FD3wg+vngg6PEcOaZ5bWaAtuiNG15vvjqU7UUlRJH1o49FpYvr7tFabvL68VX4yZSZEocWdtrryh5SKFo3ESKTOU2IiISixKHiIjEosQhIiKxKHGIiEgsGhwXSUDWZcSq6pIkKXGIJCDri7OquiRJShzS9lrVOtCnfGkXShzS9lp1Uden/OaEmGBDjClkGhwXkVSFmGBDjClkShwiIhKLEodIAdUbn8nD4pASPo1xiGQsif519ctLktTiEGmRkX7KV/+65I1aHFKTqkziS+L3Yl8dvBdLCL//0f5tZD05spYQYwpZZonDzDqAm4H9gZ3AEnf/ftU5xwO/AH5fOnSru1+eZpztSp+Cw5TW779Rchjt30bWia+WEGMKWZYtjn7gf7r7w2Y2AXjIzJa5+5qq8+5x9/dnEJ9I29IHB2kkszEOd9/o7g+Xft4CrAUOzCoeERFpThCD42Y2AzgKWFHj7neYWbeZ3WVmR6QamEgK1I8ueZP54LiZ7Qn8DLjY3as7Gh8GDnL3rWZ2MvBzYGad51kMLAaYPn16ghGLtFat/vXqQfF2oaKMfMi0xWFm44mSxi3ufmv1/e7+qrtvLf28FBhvZpNrPZe7L3H3TnfvnDJlSqJxtwNNIMtWyL//JGPT2Eo+ZFlVZcANwFp3/8c65+wPbHJ3N7P5RInuxRTDbFv6dJetrH//jcpTs45NspdlV9W7gL8DVptZV+nYF4DpAO5+LfBh4AIz6wdeB85wd88iWJEiGa5LSMlBGskscbj7vUDDjlx3vwq4Kp2IRNqHuoRkNIKoqhIRkfzIvKpKpNVUmZNfWvojH5Q4pHDUDZO+ViVrJfZ8UFeViIyaknV7UeIQaUMhzxOR8KmrSqQNqUtIRkOJQ6QBDbSLDKWuKimcVnbDqO9eZCi1OKRw1BJIn8po24sSh0jOhNh9pmTdXtRVJZIz6j6TrClxiIhILEocIg1ovoPIUBrjEGlAffetFeL4jMSnxCG5ULQLTtFeT7M0PlMM6qqSXCjaBWc0r0fdZ5I1tThEcqbILRLJByUOGaRdu1BEpHmZJg4zWwR8HxgL/NDdv1F1/67AzcDRwIvA6e7+TNpxtpOidQkViZK6hCKzxGFmY4F/At4L9AIPmtnt7r6m4rRzgZfd/RAzOwP4JnB6+tFKq+jiN3JFSOpamqQYsmxxzAfWufvTAGb2U+AUoDJxnAJ8pfTz/wOuMjNzd08zUGmdkV78inbBKdrraZY+HBRDlonjQKCn4nYvsKDeOe7eb2avAPsAL6QSoQSjaBecor0eaS9ZluNajWPVLYlmzolONFtsZivNbGVfX9+ogxMRkdqyTBy9QEfF7WnAhnrnmNk4YC/gpVpP5u5L3L3T3TunTJmSQLjtQXMERGQ4WXZVPQjMNLODgT8AZwAfqTrnduAc4D7gw8BvNb6RLHWhhKtdx0UkPJkljtKYxYXAr4jKcW9090fN7HJgpbvfDtwA/NjM1hG1NM7IKl4ZnXrVVAN08RuekrqEItN5HO6+FFhadezLFT//ETgt7bik9RolDb9MjUiRPNHMcZGc0BwYCYUWORTJiSJMAJRiUOIQEZFYlDhERCQWJQ5JheaHiBSHBsclFRq8FSkOtThEckKtNgmFWhwiOaFWm4RCiUMEzZEQiUNdVSJojoRIHGpxoE+bkhz9bUkRqcWBPm22u4lXTEzsufW3JUWkxCFtTxdxkXiUOEREJBYlDpEGNEdCZCglDpEGNIAtMpSqqtCWnFLfaKui9LclRaTEgT5VtrtGF/fRVkXpb0uKSIlDciHJ+RCNHm9ftVE9t0gRZZI4zOxbwAeAN4CngI+7++Ya5z0DbAF2AP3u3plmnBIOzYcQCUdWg+PLgCPdfQ7wBPD5Bue+x93nKWmIiIQhk8Th7ne7e3/p5v3AtCziEBGR+EIox/0EcFed+xy428weMrPFjZ7EzBab2UozW9nX19fyIKU9aQ8MkaESG+Mws18D+9e464vu/ovSOV8E+oFb6jzNu9x9g5ntCywzs8fcfXmtE919CbAEoLOz00f9AkRQVZRILYklDnc/sdH9ZnYO8H7gBHeveaF39w2l78+b2W3AfKBm4pBi03wIkXBkVVW1CLgUOM7dt9U5Zw9gjLtvKf28ELg8xTAlIPrkLxKOrMY4rgImEHU/dZnZtQBmdoCZLS2dsx9wr5l1Aw8Ad7r7L7MJV0REBmTS4nD3Q+oc3wCcXPr5aWBumnFJcWgDJZHkhFBVJdJymjAokhwlDhERiUWJQ0REYlHiEBGRWJQ4REQkFiUOKSQtFSKSHO3HIYWkkluR5KjFISIisShxiIhILEocIiISixKHiIjEosQhIiKxWJ2tMHLNzPqAZ7OOI4bJwAtZBzECijs9eYwZFHeaRhvzQe4+pZkTC5k48sbMVrp7Z9ZxxKW405PHmEFxpynNmNVVJSIisShxiIhILEocYViSdQAjpLjTk8eYQXGnKbWYNcYhIiKxqMUhIiKxKHEEwsy+ZmarzKzLzO42swOyjqkZZvYtM3usFPttZvbmrGMajpmdZmaPmtlOMwu+csbMFpnZ42a2zsw+l3U8zTCzG83seTN7JOtYmmVmHWb272a2tvT38ZmsY2qGme1mZg+YWXcp7q8m/m+qqyoMZjbR3V8t/fzfgcPd/fyMwxqWmS0Efuvu/Wb2TQB3vzTjsBoys8OAncB1wGfdfWXGIdVlZmOBJ4D3Ar3Ag8CZ7r4m08CGYWbvBrYCN7v7kVnH0wwzmwpMdfeHzWwC8BBwag5+1wbs4e5bzWw8cC/wGXe/P6l/Uy2OQAwkjZI9gFxkdHe/2937SzfvB6ZlGU8z3H2tuz+edRxNmg+sc/en3f0N4KfAKRnHNCx3Xw68lHUccbj7Rnd/uPTzFmAtcGC2UQ3PI1tLN8eXvhK9fihxBMTM/sHMeoCPAl9rVJv/AAACpUlEQVTOOp4R+ARwV9ZBFMyBQE/F7V5ycDHLOzObARwFrMg2kuaY2Vgz6wKeB5a5e6JxK3GkyMx+bWaP1Pg6BcDdv+juHcAtwIXZRls2XNylc74I9BPFnrlmYs4Jq3EsF63RvDKzPYGfARdX9QQEy913uPs8ohb/fDNLtHtQOwCmyN1PbPLUfwXuBC5LMJymDRe3mZ0DvB84wQMZNIvxuw5dL9BRcXsasCGjWAqvNEbwM+AWd78163jicvfNZvYfwCIgscIEtTgCYWYzK25+EHgsq1jiMLNFwKXAB919W9bxFNCDwEwzO9jMdgHOAG7POKZCKg0y3wCsdfd/zDqeZpnZlIFqRjN7E3AiCV8/VFUVCDP7GXAoUbXPs8D57v6HbKManpmtA3YFXiwduj/0ajAz+xvgB8AUYDPQ5e5/nW1U9ZnZycD3gLHAje7+DxmHNCwz+wlwPNGKrZuAy9z9hkyDGoaZ/SVwD7Ca6P8hwBfcfWl2UQ3PzOYANxH9fYwB/q+7X57ov6nEISIicairSkREYlHiEBGRWJQ4REQkFiUOERGJRYlDRERiUeIQSUBppdXfm9nepduTSrcPMrNfmtlmM7sj6zhFRkKJQyQB7t4DXAN8o3ToG8ASd38W+Bbwd1nFJjJaShwiyfkucIyZXQz8JfAdAHf/DbAly8BERkNrVYkkxN23m9n/An4JLCwtiy6Se2pxiCTrJGAjkIvNjESaocQhkhAzm0e0c98xwN+XdpgTyT0lDpEElFZavYZoT4f1RAPi3842KpHWUOIQScangPXuvqx0+2pglpkdZ2b3AP8GnGBmvWYW7Mq8IrVodVwREYlFLQ4REYlFiUNERGJR4hARkViUOEREJBYlDhERiUWJQ0REYlHiEBGRWJQ4REQklv8P1fNSLLEgkSoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XuUXGWZ7/Hvk3TCLQmEJEBCOgalBRQEhjZB5RIWiklEYWa8BG8c9UxgljB65jiIehTFmYPGy+iBQYjCEWcQZzyKsjQoUc8scB0JBAyKAumQW4dEEhLIxQRIp5/zx66mq6qrqmt319773bt+n7VqVdWuXdVPdXXvp973ed93m7sjIiLSrDFZByAiIvmixCEiIrEocYiISCxKHCIiEosSh4iIxKLEISIisShxiIhILEocIiISixKHiIjE0pF1AEmYOnWqz549O+swRERy46GHHnrG3ac1s28hE8fs2bNZuXJl1mGIiOSGmW1odl91VYmISCxKHCIiEosSh4iIxJJ44jCzW81sq5k9Wrbts2b2lJmtKl0W1nnufDN7wszWmNnVSccqIiLDS6PF8W1gfo3t/+zup5Uuy6ofNLOxwL8AC4BXAZeY2asSjVRERIaVeOJw93uBHSN46hxgjbuvdfcXge8BF7U0OBERiS3LGscVZva7UlfW5BqPHwv0lt3fVNomIiIZyipxfAN4BXAasAX4So19rMa2uue5NbPFZrbSzFZu27atNVGK5FhvL1x5JcyZE1339g7/HJFmZDIB0N2fHrhtZt8EflJjt01AZ9n9mcDmBq+5FFgK0N3drROpS1vr7YVTT4U9e2D/fli1Cm6/HR55BDo7h3++SCOZtDjMbHrZ3b8EHq2x24NAl5kdZ2bjgUXAXWnEJ5J3S5YMJg2IrvfsibaLjFbiLQ4zuwOYB0w1s03ANcA8MzuNqOtpPXBZad8ZwLfcfaG795nZFcDPgbHAre7+h6TjFSmCFSsGk8aA/fvhgQeyiUeKJfHE4e6X1Nh8S519NwMLy+4vA4YM1RWRxubOjbqnypPHuHFRvUNktDRzXKSArroKJkyIkgVE1xMmRNtFRkuJQ6SAOjujQvhll0WtjMsuU2FcWqeQy6qLFEVvb1TQXrEi6n666qrmD/6dnXD99cnGJ+1JiUMkUBpSK6FSV5VIoDSkVkKlxCESKA2plVApcYgEpHyZkOefh46qzuRQhtRqOZP2phqHSCCqaxodHXDgQHTd1xfOkFrVXkQtDpFAVNc0+vqipHHSSWENqVXtRdTiEAlEvZrGIYdEj4VCtRdRi0MkEHPnDs70HhBKTaNcXuKU5ChxiAQiL8uE5CVOSY4Sh0gg8rJMSBpxatRW2My9eOc86u7u9pUrV2YdhoiMQPWorYEWTYhJtEjM7CF3725mX7U4RCQoGrUVPiUOEQmKRm2FT4lDRIKiUVvhU+IQkZqyKlAnOWpLRffWUHFcRIbIukA9cB6S++6D/n4YMwbOPjve+UhqvaaK7vUFVRw3s1vNbKuZPVq27Utm9riZ/c7M7jSzI+o8d72Z/d7MVpmZMoFISrIuUHd2Rkli40Z4/PHo4H7zzdGBf6SthKzfU5Gk0VX1bWB+1bblwMnu/hpgNfCJBs8/z91PazYTikilkXTPhFCgbvWBPoT3VBSJJw53vxfYUbXtHnfvK929H5iZdBwi7Wige+bmm+HBB5v/1h5CgbrVB/oQ3lNRhFAc/yBwd53HHLjHzB4ys8WNXsTMFpvZSjNbuW3btpYHKZJHn/40PPdc/G/tISwr0uoDfQjvqSgyTRxm9imgD7i9zi5vcPe/ABYAHzazc+q9lrsvdfdud++eNm1aAtFK0jTipbV6e+Hf/g2qx7808609hOVPWn2gD+E9FUUqo6rMbDbwE3c/uWzbpcDlwPnuvreJ1/gssMfdvzzcvhpVlT8a8dJ6V14JN94YjUoqZwbvfz9MnBh1B82dO7rRSkkaGF31wAPRwT7UOIsgzqiqTM7HYWbzgY8D59ZLGmZ2GDDG3XeXbl8AXJtimJKiRoXQ66/PNra8WrFiaNKAKHH86Eewd2/4Z/Dr7NTnH6I0huPeAfwGOMHMNpnZh4AbgInA8tJQ25tK+84ws2Wlpx4N/NrMHgEeAH7q7j9LOl7Jhka8tF6tGsGYMdDVNZg0QMNSJb7EWxzufkmNzbfU2XczsLB0ey1waoKhSUDmzo2++ZYnj5BHvAx0oYTc1XPVVVFLorr77+CDlaRldEIYVSWSqxEvIx3imrZ6xeCzz9awVBkdLTkiwchLIfTKK6NkUd06uuyyfPTHhzQQobrl9t73RiPBQm7JFVWc4rgSh7Sd3t5ofsPdpdlDCxbA5z/f/AFqzpyopVFr+4oVrYszSSEk6eoE1tEBBw7A2LHQ16eRdWkLflSVSFZ6e+GUU2DnzsFtt90WjTL6/e+bO0DlrR5TSwijlapH0vX1VV5rZF24VOOQQmh28uCSJbBr19Dtu3c3P6ooT/WYkNUaSVdNRfswqcUhuVfd5dFoXsKKFUNnUkM036HZA9RA0Tnrrp68q9Vyq5a3lly7UItDci/OKqpz50YT4KqNGRPvADXQ1bNiRXStpBFfdcutoyP6bDpKX2fVkguXEofkXpzJg1ddBZMmDd0+cWJ7HKBCWg+serjw5ZfDb34TXWstqbCpq0pyL06xurMzKoI3GlUV8uS+0cRW3aX329/CN78JJ544+rPrjVStIv3cuenGIPFpOK7kXivnJYQ0x6HVsdWafzIgpPcp2Qjq1LEiSWvlctkhn150tLE1GsUU0vuU8ClxSCG0qlidxGKLraorjDa2WosexnmtkOojki3VOETKtHpyX5yhwknHVr3oYbVGr9XK9yH5pxaH5FqrvwW3enJfK7u+RhtbeZfeqafCQQc1P/Q15C48SZ8Sh+RWEqvUtvr0oqPpXqpOijD62Aa69Fatgp6e5oe+6nwpUk5dVZJbSZ01sNE6TnGHw77qVfDQQ5Vn4mume6lR11Cr1m2Ks15VEdbnktZRi0NyK+1vwXFbOL290eKJ1advPfTQ4buXQusaKsr6XCrwt4YSh+RWrVFCSX4LjnswX7IkOkVrOTO4+OLhu5dC6xpqdRdeFvJyAq48SCVxmNmtZrbVzB4t23akmS03s57S9eQ6z720tE+PmV2aRrySD2l/C457MK+1vzs89tjwPyvtpNiMvK/PFVorLs/SanF8G5hfte1q4Jfu3gX8snS/gpkdCVwDzAXmANfUSzDSftL+Fhz3YD6ag39RuoZCElorLs9SSRzufi+wo2rzRcBtpdu3ARfXeOqbgeXuvsPdnwWWMzQBSRtL81tw3IP5aA7+aSXFdurzD7EVl1eprVVlZrOBn7j7yaX7z7n7EWWPP+vuk6ue8zHgYHf/x9L9TwP73P3LjX6W1qqSpMQ95WoIp2itJ+R1uZLQbu83riKdOrbGmROomenMbDGwGGDWrFlJxiRtLO4pV9M4RetIV8xNajhzqHQCrtbJMnE8bWbT3X2LmU0HttbYZxMwr+z+TOA/a72Yuy8FlkLU4mhtqCJhGs1SICH0+TeT9Fq5zH0I51ovgiyH494FDIySuhT4cY19fg5cYGaTS0XxC0rbJDDt1Fceknqthre8ZfjPot6ih/v2pfP5NTM8VkNoA+XuiV+AO4AtwH6iVsSHgClEo6l6StdHlvbtBr5V9twPAmtKlw808/POOOMMl/Rs3Og+ebL7uHHuEF1Pnhxtl2S99rXR77zepdFnMfC5dXRUPqejI53P74orBv9myuO94op4+0hrACu9yWN6WqOqLnH36e4+zt1nuvst7r7d3c93967S9Y7Svivd/b+WPfdWdz++dPnfacQr8Wh8fHaaWSq91mcx0P3T2Tn0VLp9fel8fs10ld17b/bdaTJU6MVxyYEQ+srb1XBLpcPQz6K6LtLMc5Iw3PpXvb3wxBNDn9fRoSG0WdOSIzJqGh+fner5HqecMrhU+oDqz6K6hVhLGp/fcPNcliyBAweGPm/sWE2EzJoSh4xavQPAe9+b74J5CAX/ZmIonwT505/CxImNJx02OoVsveckYbhJjitWRN1m1U48UUNoM9dsMSRPFxXH07dxY1SwnDMnur7//nwXzEMo+I80hurPonr/egXnU06p/5wsqDCeLmIUxzM/yCdxUeLIXt7/6UOIP6kYQkiKzchLnEURJ3Goq0oSkfeCeQjxJxVDXpZIz0uc7UijqiQReT9jXAjxNxvDSGZW52UGdV7ibDepLXKYJi1ymL28LygXQvzNxBBCnFIMcRY5VFeVJCLv3QwhxN9MDJp82Wb27Imaod//fvTHkBG1OERybM6caA2nWttXrEg/HmmBffvgySehp2foZfPmwf2uvhquu45J101i94u7h7zMxPET2fWJXU3/2CItqy4iDbSiFtOqA4/E8OKLsHZtZVJYvTq6rp6sM20adHXBBRdE1+UXqPnZNdreCkocIjlWveTISCbvZXHgaQv798P69bVbDhs2QH//4L5HHhklgnPPHZocDj88s7dQjxKHSI7p5EQZO3AANm4cbC2UX9avr5z6PmlSlAjmzo2WVShPDlOmZPYWRkKJQyQBrTz50HA0ZDVh/f3w1FOV3UkDl7Vro26nAYcdFiWC00+Hd76zMjlMmwZW66Sm+aPEIdJiozkrn2TEHbZsqd2ttGYNPP/84L4HHwzHHw8nnQRve1uUFF75yuj6mGMKkxwaUeIQabF2O5d3brjDtm21Ww5r1sCf/zy47/jx8PKXVxalB5LDscfCmHBmMkwcP7Hu4IakKHGItFgIy5XEkcWBJ1Hbt9duOfT0wK6yUWIdHXDccVEymDevsltp1qxo/fYcyGLkmxKHSIuFsFxJHLkccrtzZ+2CdE8PPPvs4H5jxrDucGf1kU7PCdAzBVZPgZ4jYcdRB7Pjf6zO7j3kWGaJw8xOAP69bNPLgc+4+9fK9pkH/BhYV9r0Q3e/NrUgRUagFUNki2DUAwT27Knfcti2bXA/s+iFu7rgXe+qbDkcdxwv/8LBtV//wJ5Rvb92llnicPcngNMAzGws8BRwZ41d73P3C9OMTWQ0NEQ2xgCBvXtrz5JevRr+9KfKF50xI0oGF19cmRxe8YqoYC2pCaWr6nzgSXffkHUgIq3Q7kNkywcIjOcFXrH/SU7a1cPKRT10nlyWIDZtqnzi0UdHyWDBgqHJYcKEbN6MDBFK4lgE3FHnsdeZ2SPAZuBj7v6H9MISkabs3w/r1r2UEM75jx7eur+HLnqYxUbG0g8HgP8HPDElSgbnnTd0lvSkSVm/E2lC5onDzMYDbwM+UePhh4GXufseM1sI/AjoqvM6i4HFALNmzUooWpE21tcXLZVRq+awfn00i7rkLeMP54/WxW/8ddzGpfTQxbqxXZz1gS6WfHNydu9BWiLzxAEsAB5296erH3D3XWW3l5nZjWY21d2fqbHvUmApRKvjJhmwZCfNGdltqb8/+iXXWnxv3brKoWITJkSthDPOgEWLKloO2/dN5YLTbMgAgX//TPpvqXDDjQMQQuK4hDrdVGZ2DPC0u7uZzSE6f8j2NIOTcGhG9shVrIDrMGM3dG2HU3YdxPXH/91ggnjySXjhhcEnHnJIlAxOOQX+6q8qu5WOPrruLOlOwhkgkMvhxoHLNHGY2aHAm4DLyrZdDuDuNwFvB/7WzPqAfcAiL+IJRKQpmpEdgzs8/fRLrYZPLNtN13bo2gHH74DDXmo4vAAH/a+o+NzVBQsXViaHGTNqzpJuZin2dh8gUGSZJg533wtMqdp2U9ntG4Ab0o6rKIp2noW8zchOQnXLYcreKBm8ZufB3HzixwZbDmvWwO7Bz/5jY2Dt5Gji26+Oi64HJsNt+PKfY8+S1lLs7S2EripJSNH+ufM2I7tlnn32pZbD3/98N107eKn1MPmltfeehzH/c3AJjbPOqmg5HHLb8RyolxtysrSGhEOJQ3Kj0DOyd+8eWoweuGwfLOt9Bth4eNRauOPk6Hqg9fDEV/ZFi/PVUDdpiIyAEofkRu5nZP/5z1EXUq3hrE9XDSqcOTNqLfz1X1e0HA797qt5YVyd16+TNERaTYlDciWtguuIh/0+/3ztJTR6eqKTAZU75pgoIbzlLYPJ4ZWvjArVhx5a8+XrJo1haEhq9opUc1TiEKky7LDfF1+smCVd0bXU2xuNaBowdWqUEN74xsrRSscfDxPTO2i3+sCkRBRfkWqOShwFpn/ukVmyBPbt7mNW33q66KFrfw8n7Oxh91k90FGaJd3fP/iEyZOjZHD22UOX0DjiiJbGFspnmrdvyNJaShwFpn/uYRw4ELUQqorR//DLHr7at45x9L20667+iWx+pgve+lp497sru5amTGnwQ1pLn6mEQIlDiq2/P6ot1Ko5PPlk1O004NBDo+UyZp7KHevezuP9XUTL9HXxbMdRXPZBS21CWxH6w4vwHqQ2JQ4JxogL0u7RuRtqJYc1a2DfvsF9Dzooqi+ccAJceOHgeaS7umD6dDBjai988dRsh/0WoT+8CO9BamuYOMxsEjDN3Z+s2v4ad/9dopFJWxm2IO0OzzxT+3Sha9ZETxwwbtzgEhpvelNlzWHmzJpLaJTL47BffbsPXyj1qVaomzjM7J3A14CtZjYO+C/u/mDp4W8Df5F8eNIuBtahmrB/x0sF6RN39vD0+T10Hl5KEDt3Dj5h7NjBWdLnnFNZc5g1a9SzofO2zpK+3YevSAm8UYvjk8AZ7r6ltDLtv5rZJ939h0DtJTFFmrFz55BWw+I7e/js/h6msOOl3fr7jT9teBmc2wXveU9lt9Ls2VHLIiP6ht/YpOt0QqYia5Q4Otx9C4C7P2Bm5wE/MbOZgFaolcb27Kk9S3r1ati2rXLfzk7GT+7iB/vewRNlBenejpfzgcUHBfnNX9/wG9PvodgaJY5dZvaKgfpGqeUxj+gsfK9OIzgJ3L599ZfQ2LKlct/p06MWw0UXDT2X9CGHcGgvXJ1xQTokReoPr1aE99DuGiWOj1PVJeXuu81sPrVP8ypF9MILsHZt7eTQ21u571FHRcngzW8eOkt6woSGPyaPBekkFbm7K633pu7E5DRKHLcBN5vZV9y9D8DMjga+ApwAXJtCfJKG/fuj2dC1ksOGDZWzpKdMiRLBuecOFqMHksPhh48qjLwVpENS5BbKSKk7MTmNEscZwHXAb83sI8ApwN8DS4D3pxCbtNKBA1ESqJUc1q+HvsFZ0kyaFCWDM8+E972vsvVw5JGZvQWpT9+gJU11E4e7PwtcXkoavwA2A2e6+6a0gpOY+vth06ahxeienqi7qfwMSIcdFiWC00+Hd76zMjlMm1b3XNIS0Tf8xvT7KbZG8ziOAL4IzAXmAwuBu83sI+7+q5Tik2rusHlz/SU0nn9+cN9DDom6kF79arj44srkcMwxSg6joG/4jen3U2yNuqoeBm4EPlyqcdxjZqcBN5rZBne/pBUBmNl6YDdwAOhz9+6qxw34OlHi2ks0EfHhVvzsYLnD1q21k0NPD+zdO7jv+PGDs6QXLKhMDjNmDDtLWiRP0ih4q6g+vEaJ45zqbil3XwW83sz+psVxnOfuz9R5bAHQVbrMBb5Rui6GLVtg+fKhyWF32R9uR8fgLOl58yqTQwtmSYukbaQH5zgF75F2l6moPrxGNY66tQx3/2Yy4dR0EfAdd3fgfjM7wsymD0xOzL1HHoFLL41aBrNnR8ng9a+vTA6zZ0fJQ6Qg0jg4q3WQnBCORk7UDebAze6+tOrxY4HyCQObStsqEoeZLQYWA8yaNSu5aFvtrLPg8cejFoXOGV0o6vKQogohcbzB3Teb2VHAcjN73N3vLXu8VgV3yJInpYSzFKC7uzs/S6JMmBAt8S015fngqy4PKarMK6fuvrl0vRW4E5hTtcsmoHz+8EyiocHSBnTwFQlPponDzA4zs4kDt4ELgEerdrsLeL9FzgR2Fqa+ISKx1Ctst3J+SBo/I++y7qo6GrgzGnFLB/Bdd/+ZmV0O4O43AcuIhuKuIRqO+4GMYpWCynN3WF6NdMRTGp+HPvPhZZo43H0tcGqN7TeV3Xbgw2nGJe1F3WHp08E53zKvcYgUlbo8pKiy7qoSaSjPax7pW7UUlRKHBC2vB1/VTaTI1FUlkgDVTaTIlDhERCQWJQ4REYlFiUNERGJR4hARkVg0qkokAVkPI9aoLkmSEodIArI+OGtUlyRJiUPaXqtaB/qWL+1CiUPaXqsO6vqW35wQE2yIMYVMxXERSVWICTbEmEKmxCEiIrEocYgUkFbmlSSpxiGSsST619UvL0lSi0OkRUb6LV/965I3anFITRplEl8Svxf7nFXcD+H3P9q/jawnR9YSYkwhyyxxmFkn8B3gGKAfWOruX6/aZx7wY2BdadMP3f3aNONsV/oWHKa0fv+NksNo/zayTny1hBhTyLJscfQB/93dHzazicBDZrbc3f9Ytd997n5hBvGJtC19cZBGMqtxuPsWd3+4dHs38BhwbFbxiIhIc4IojpvZbOB0YEWNh19nZo+Y2d1m9upUAxNJgfrRJW8yL46b2QTgB8BH3b26o/Fh4GXuvsfMFgI/ArrqvM5iYDHArFmzEoxYpLVq9a9XF8XbhQZl5EOmLQ4zG0eUNG539x9WP+7uu9x9T+n2MmCcmU2t9VruvtTdu929e9q0aYnG3Q40gSxbIf/+k4xNtZV8yHJUlQG3AI+5+1fr7HMM8LS7u5nNIUp021MMs23p2122sv79NxqemnVskr0su6reALwP+L2ZrSpt+yQwC8DdbwLeDvytmfUB+4BF7u5ZBCtSJMN1CSk5SCOZJQ53/zXQsCPX3W8AbkgnIpH2oS4hGY0gRlWJiEh+ZD6qSqTVNDInv7T0Rz4ocUjhqBsmfa1K1krs+aCuKhEZNSXr9qLEIdKGQp4nIuFTV5VIG1KXkIyGEodIAyq0iwylriopnFZ2w6jvXmQotTikcNQSSJ+G0bYXJQ6RnAmx+0zJur2oq0okZ9R9JllT4hARkViUOEQa0HwHkaFU4xBpQH33rRVifUbiU+KQXCjaAado76dZqs8Ug7qqJBeKdsAZzftR95lkTS0OkZwpcotE8kGJQyq0axeKiDQv08RhZvOBrwNjgW+5+xeqHj8I+A5wBrAdeJe7r087znZStC6hIlFSl1BkljjMbCzwL8CbgE3Ag2Z2l7v/sWy3DwHPuvvxZrYI+CLwrvSjlVbRwW/kipDUtTRJMWTZ4pgDrHH3tQBm9j3gIqA8cVwEfLZ0+/8AN5iZubunGai0zkgPfkU74BTt/TRLXw6KIcvEcSzQW3Z/EzC33j7u3mdmO4EpwDOpRCjBKNoBp2jvR9pLlsNxrca26pZEM/tEO5otNrOVZrZy27Ztow5ORERqyzJxbAI6y+7PBDbX28fMOoDDgR21Xszdl7p7t7t3T5s2LYFw24PmCIjIcLLsqnoQ6DKz44CngEXAu6v2uQu4FPgN8HbgV6pvJEtdKOFq17qIhCezxFGqWVwB/JxoOO6t7v4HM7sWWOnudwG3AP9qZmuIWhqLsopXRqfeaKoBOvgNT0ldQpHpPA53XwYsq9r2mbLbzwPvSDsuab1GScOvUSNSJE80c1wkJzQHRkKhRQ5FcqIIEwClGJQ4REQkFiUOERGJRYlDUqH5ISLFoeK4pELFW5HiUItDJCfUapNQqMUhkhNqtUkolDhE0BwJkTjUVSWC5kiIxKEWB/q2KcnR35YUkVoc6Ntmu5t03aTEXlt/W1JEShzS9nQQF4lHiUNERGJR4hBpQHMkRIZS4hBpQAVskaE0qgqdklPqG+2oKP1tSREpcaBvle2u0cF9tKOi9LclRaTEIbmQ5HyIRs+3z9moXlukiDJJHGb2JeCtwIvAk8AH3P25GvutB3YDB4A+d+9OM04Jh+ZDiIQjq+L4cuBkd38NsBr4RIN9z3P305Q0RETCkEnicPd73L2vdPd+YGYWcYiISHwhDMf9IHB3ncccuMfMHjKzxY1exMwWm9lKM1u5bdu2lgcp7UnnwBAZKrEah5n9AjimxkOfcvcfl/b5FNAH3F7nZd7g7pvN7ChguZk97u731trR3ZcCSwG6u7t91G9ABI2KEqklscTh7m9s9LiZXQpcCJzv7jUP9O6+uXS91czuBOYANROHFJvmQ4iEI6tRVfOBjwPnuvveOvscBoxx992l2xcA16YYpgRE3/xFwpFVjeMGYCJR99MqM7sJwMxmmNmy0j5HA782s0eAB4CfuvvPsglXREQGZNLicPfj62zfDCws3V4LnJpmXFIcOoGSSHJCGFUl0nKaMCiSHCUOERGJRYlDRERiUeIQEZFYlDhERCQWJQ4pJC0VIpIcnY9DCklDbkWSoxaHiIjEosQhIiKxKHGIiEgsShwiIhKLEoeIiMRidU6FkWtmtg3YkHUcMUwFnsk6iBFQ3OnJY8yguNM02phf5u7TmtmxkIkjb8xspbt3Zx1HXIo7PXmMGRR3mtKMWV1VIiISixKHiIjEosQRhqVZBzBCijs9eYwZFHeaUotZNQ4REYlFLQ4REYlFiSMQZvZ5M/udma0ys3vMbEbWMTXDzL5kZo+XYr/TzI7IOqbhmNk7zOwPZtZvZsGPnDGz+Wb2hJmtMbOrs46nGWZ2q5ltNbNHs46lWWbWaWb/18weK/19fCTrmJphZgeb2QNm9kgp7s8l/jPVVRUGM5vk7rtKt/8OeJW7X55xWMMyswuAX7l7n5l9EcDdP55xWA2Z2UlAP3Az8DF3X5lxSHWZ2VhgNfAmYBPwIHCJu/8x08CGYWbnAHuA77j7yVnH0wwzmw5Md/eHzWwi8BBwcQ5+1wYc5u57zGwc8GvgI+5+f1I/Uy2OQAwkjZLDgFxkdHe/x937SnfvB2ZmGU8z3P0xd38i6ziaNAdY4+5r3f1F4HvARRnHNCx3vxfYkXUccbj7Fnd/uHR7N/AYcGy2UQ3PI3tKd8eVLokeP5Q4AmJm/2RmvcB7gM9kHc8IfBC4O+sgCuZYoLfs/iZycDDLOzObDZwOrMg2kuaY2VgzWwVsBZa7e6JxK3GkyMx+YWaP1rhcBODun3L3TuB24Ipsox00XNylfT4F9BHFnrlmYs7Mx4WOAAACT0lEQVQJq7EtF63RvDKzCcAPgI9W9QQEy90PuPtpRC3+OWaWaPegzgCYInd/Y5O7fhf4KXBNguE0bbi4zexS4ELgfA+kaBbjdx26TUBn2f2ZwOaMYim8Uo3gB8Dt7v7DrOOJy92fM7P/BOYDiQ1MUIsjEGbWVXb3bcDjWcUSh5nNBz4OvM3d92YdTwE9CHSZ2XFmNh5YBNyVcUyFVCoy3wI85u5fzTqeZpnZtIHRjGZ2CPBGEj5+aFRVIMzsB8AJRKN9NgCXu/tT2UY1PDNbAxwEbC9tuj/00WBm9pfA9cA04Dlglbu/Oduo6jOzhcDXgLHAre7+TxmHNCwzuwOYR7Ri69PANe5+S6ZBDcPMzgLuA35P9H8I8El3X5ZdVMMzs9cAtxH9fYwB/sPdr030ZypxiIhIHOqqEhGRWJQ4REQkFiUOERGJRYlDRERiUeIQEZFYlDhEElBaaXWdmR1Zuj+5dP9lZvYzM3vOzH6SdZwiI6HEIZIAd+8FvgF8obTpC8BSd98AfAl4X1axiYyWEodIcv4ZONPMPgqcBXwFwN1/CezOMjCR0dBaVSIJcff9ZvYPwM+AC0rLoovknlocIslaAGwBcnEyI5FmKHGIJMTMTiM6c9+ZwH8rnWFOJPeUOEQSUFpp9RtE53TYSFQQ/3K2UYm0hhKHSDL+Btjo7stL928ETjSzc83sPuD7wPlmtsnMgl2ZV6QWrY4rIiKxqMUhIiKxKHGIiEgsShwiIhKLEoeIiMSixCEiIrEocYiISCxKHCIiEosSh4iIxPL/ARDC/QIHKQlvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dataMat,labelMat = loadDatSet('C:/Users/exquisite/Downloads/machinelearninginaction-master/Ch05/testSet.txt')\n",
    "weights_1 = gradientDescent(dataMat,labelMat)\n",
    "weights_2 = stocGradDescent(dataMat,labelMat)\n",
    "plotBestFit(weights_1,dataMat,labelMat)\n",
    "plotBestFit(weights_2,dataMat,labelMat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "分类器"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def classifyVector(inX,weights):\n",
    "    prob = sigmoid(sum(inX*weights))\n",
    "    if prob > 0.5:\n",
    "        return 1.0\n",
    "    else:\n",
    "        return 0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "测试函数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def colicTest():\n",
    "    frTrain = open('C:/Users/exquisite/Downloads/machinelearninginaction-master/Ch05/horseColicTraining.txt') # 载入训练集\n",
    "    frTest = open('C:/Users/exquisite/Downloads/machinelearninginaction-master/Ch05/horseColicTest.txt') # 载入测试集\n",
    "    trainingSet = [] # 用于存放训练数据\n",
    "    trainingLabels = [] # 用于存放训练标签\n",
    "    # 整理数据\n",
    "    for line in frTrain.readlines():\n",
    "        currLine = line.strip().split('\\t')\n",
    "        lineArr = []\n",
    "        for i in range(21):\n",
    "            lineArr.append(float(currLine[i]))\n",
    "        trainingSet.append(lineArr)\n",
    "        trainingLabels.append(float(currLine[21]))\n",
    "    trainWeights = stocGradDescent(trainingSet,trainingLabels,500) # 得到回归系数\n",
    "    errorCount = 0 # 错误数\n",
    "    numTestVec = 0.0 # 存放测试数据个数\n",
    "    for line in frTest.readlines():\n",
    "        numTestVec += 1.0\n",
    "        currLine = line.strip().split('\\t')\n",
    "        lineArr = []\n",
    "        for i in range(21):\n",
    "            lineArr.append(float(currLine[i]))\n",
    "        if int(classifyVector(lineArr,trainWeights)) != int(currLine[21]): # 分类错误\n",
    "            errorCount += 1\n",
    "    errorRate = (float(errorCount)/numTestVec) # 计算错误率\n",
    "    print('the error rate of this test is: %f' % errorRate)\n",
    "    return errorRate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "多次测试函数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def multiTest():\n",
    "    numTests = 10\n",
    "    errorSum = 0.0\n",
    "    # 计算十次取平均\n",
    "    for k in range(numTests):\n",
    "        errorSum += colicTest()\n",
    "    print('after %d iterations the average error rate is: %f' % (numTests,errorSum/float(numTests)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the error rate of this test is: 0.388060\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.3880597014925373"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "colicTest()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the error rate of this test is: 0.358209\n",
      "the error rate of this test is: 0.358209\n",
      "the error rate of this test is: 0.283582\n",
      "the error rate of this test is: 0.432836\n",
      "the error rate of this test is: 0.253731\n",
      "the error rate of this test is: 0.417910\n",
      "the error rate of this test is: 0.253731\n",
      "the error rate of this test is: 0.477612\n",
      "the error rate of this test is: 0.552239\n",
      "the error rate of this test is: 0.507463\n",
      "after 10 iterations the average error rate is: 0.389552\n"
     ]
    }
   ],
   "source": [
    "multiTest()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
